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We show theoretically that periodic density patterns are stabilized in two counter-propagating 
Bose-Einstein condensates of atoms in different hyperfine states under Rabi coupling. In the presence 
of coupling, the relative velocity between two components is localized around density depressions in 
quasi-one-dimensional systems. When the relative velocity is sufficiently small, the periodic pattern 
reduces to a periodic array of topological solitons as kinks of relative phase. According to our 
variational and numerical analyses, the soliton solution is well characterized by the soliton width 
and density depression. We demonstrate the dependence of the depression and width on the Rabi 
frequency and the coupling constant of inter-component density-density interactions. The periodic 
pattern of the relative phase transforms continuously from a soliton array to a sinusoidal pattern 
as the period becomes smaller than the soliton width. These patterns become unstable when the 
localized relative velocity exceeds a critical value. The stability-phase diagram of this system is 
evaluated with a stability analysis of countersuperflow, by taking into account the finite-size-effect 
owing to the density depression. 

PACS numbers: 03.75.Lm, 67.85.Be, 03.75.Mn, 05.45.Yv 


I. INTRODUCTION 

In superfluid systems, counterflow of two interpene¬ 
trating fluid components is stable in the presence of a 
frictionless superfluid component. Such a flow state was 
hrst realized after the discovery of He-II, which consists of 
normal-fluid and superfluid components in the two-fluid 
model [I|. In this system, a temperature gradient causes 
a counterflow of two components along the gradient, 
termed thermal counterflow 0. Thermal counterflow 
is an important system in the field of low-temperature 
physics tha^ertain to the visualization of quantum tur¬ 
bulence IJ . Another interesting example is the coun¬ 
terflow of two superfluid components, called countersu¬ 
perflow. Although countersuperflow itself must be a fun¬ 
damental flow state in multi-component superfluid sys¬ 
tems, its study received less attention until a recent ex¬ 
perimental study on the instability of countersuperflow 
by Hammer et al. 

In their experiment it was found that novel soli¬ 
ton dynamics in quasi-one-dimension arise from the coun- 
tersuperflow instability (CSI) in miscible two-component 
BECs In this experiment, relative motion was in¬ 

duced by utilizing the Zeeman energy shift between the 
two components under a magnetic-field gradient. When 
relative velocity exceeds a critical value, countersuper¬ 
flow becomes dynamically unstable, creating soliton in 
quasi-one dimensional systems. Very recently, soliton dy¬ 
namics in a similar system 0 with an internal Joseph- 
son effect ac ting between two components, called Rabi 
coupling (ll| - ll4j | was observed. In the experiment 0, 
the magnetic-field gradient induced a spatial dependency 
of the detuning of the coupling, which played a domi¬ 
nant role in soliton nucleation. This system is interesting 
in the sense that a countersuperflow system under Rabi 
coupling can be realized if the relative motion between 
condensates is realized in a similar manner as the earlier 


experiment @ . It is expected that different kinds of soli¬ 
tons and instability developments appear in systems of 
Rabi-coupled countersuperflow. 

Motivated by these experiments, we study theoreti¬ 
cally Rabi-coupled countersuperflow of miscible binary 
condensates in quasi-one-dimensional systems. We found 
that the soliton patterns, which are distinctly different to 
those observed in Ref. [1, 0, are stabilized in this sys¬ 
tem. In the specific limits presented here, the soliton re¬ 
duces to a domain wall of a relative phase, as predicted by 
D. T. Son et al. [l^. This kind of structure is known to 
be stabilized between two vortices in the vortex-molecule 
structure in rotating Rabi-coupled two-component BECs 
However, the stability of such a structure has 
never been explored quantitatively, even for the quasi- 
one dimensional system. In this work, we present spatial 
profiles and stability-phase diagrams of the soliton pat¬ 
terns by varying the Rabi frequency, the inter-component 
coupling constant, and the relative velocity between the 
two components. 

This paper is organized as follows. Section [IT] is de¬ 
voted to the introduction of the stability analysis of Rabi 
coupled BECs and countersuperflows. In Sec. Ell vari¬ 
ational and numerical analyses are performed for the 
single-soliton solution in the limit of small relative ve¬ 
locities. Then we present the stability-phase diagram of 
the single-soliton. In Sec. lYl the problem is generalized 
to the case of multi-soliton solutions with larger relative 
velocities. Finally, in Sec. El our results are summarized 
and additional discussions are made. 


II. BASIC STABILITY ANALYSES 

Before discussing the soliton solution, we have to in¬ 
troduce bulk state, which is realized in bulk far from the 
soliton. First, we will formulate the stability of Rabi- 
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coupled two-component condensate BECs without a rel¬ 
ative velocity. Then, we will present the stability of a 
countersuperflow without Rabi coupling. 


A. Stability of Rabi-coupled condensates 


Binary BECs at zero temperature are described by 
the condensate wave function i/'i = 

(j = 1,2) in the Gross-Pitaevskii (GP) model [3- The 
Lagrangian of this system under Rabi-coupling is written 
as 


L= I -E (1) 


with the energy functional 




i=i 
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di’j 


dx 










( 2 ) 


where mj is the atomic mass, and jij is the chemical 
potential of the jth component. The coefficient gjk = 
2'K’h?ajk/mjk of the density-density interaction is repre¬ 
sented by the effective mass mjk = (^ 7 ^ + ^^.d 

the s-wave scattering length ajk = Ofcj between the jth. 
and Ath components. The last term on the ri ght hand 
side of Eq. (I2|) represents the Rabi coupling [ll|. We 
may set the Rabi frequency 17 as 17 > 0 without loss of 
generality. 

In the presence of Rabi coupling, population transfer 
occurs between two components, and as such we set m-i = 
m 2 = m and /ii = /i 2 = We restrict ourselves to 
the case of gu = 322 = g > 0 and 512 > 0 , which is 
typically satisfied in Rabi-coupled condensates 
Without Rabi coupling, miscible states of binary BECs 
are unstable for g < 512 , and the condensates undergo 
a phase separation. Hence, we use the non-dimensional 
variable 



(3) 


as a characteristic parameter of this system. 

The ground state is obtained by neglecting the spatial 
dependence of the order parameters i/'j ■ As such, Eq. ([5]) 
reduces to 


E = y |e + E \djknjnk 




— hQ.^nin 2 cos 


(4) 


where the relative phase is 


e_ = 61 - 62 . (5) 

The ground state is obtained by minimizing Eq. ([5]) with 
respect to 6-, ni, and 77 , 2 . By using 


and 


9(1 + 1) 



17o = 


2 gno 

h 


the ground state with 0 _ = 0 is written as 


ni = 772 = no 

for 7 < 1 -|- 17/17o and as 


( 6 ) 

(7) 

( 8 ) 


ni,2 = no 



/y - 1 - ll/llo\ 

7^1 ) 


(9) 


for 7 > 1 -I- 17/17o [i3| • The state ([5]) is the bulk state of 
our soliton solution, by assuming the condition 


17 ^ 17° 
17o 17o 


( 10 ) 


For 17 = 0, the criterion (nni) reduces to that of the phase 
separation; 7 = 1. This criterion represents the global 
energetic stability since the analysis is based on the com¬ 
parison between the energies of states dH) and ®. 

The local stability, the so-called linear stability, of 
the Rabi-coupled condensates was inv estig ated using the 
Bogoliubov-de Gennes (BdG) theory Here, we 

investigate the linear stability around the bulk state ([ 8 ]) ; 
'4’j — \/no- By linearizing the equation of motion ob¬ 
tained from Eq. © with respect to a collective per¬ 
turbation 5'tpj{x,t) = tpj[x,t) — y/nQ = — 

and diagonalizing the linearized equations, 
we obtain the dispersion relations 


q^e fq^e 


17f 


2 V 2 


-I-1 -I- 7 


( 11 ) 


and 
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17o 


q^e . n \ fq^^' 
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where we used 




y/2mgno 


(13) 


Since the system is unstable for < 0, we find from 
Eq. (fT^ that the stability condition of state ([ 8 ]) is repre¬ 
sented again by Eq. (fTUl) . In this work, we consider the 
parameter region that satisfies condition m, as shown 
in Fig. [11(a). 
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(a) (b) 




FIG. 1. The stability-phase diagram of the bulk state in 
Rabi coupled condensates (a) and a uniform countersuperflow 
(b). The phase boundaries in (a) and (b) represent Eqs. (IIOII 
and GSl), respectively. 


B. Stability of the countersuperflow 


Before we discuss the stability of the Rabi-coupled 
countersuperflow, it is useful to demonstrate the linear 
stability of the countersuperflow. A stationary solution 
of our system is described by the time-independent GP 
equation 


f -h^ 9^ I , l2 I , |2\ / / 

I + 93km\ J ipj - -Y^k 

(14) 


for k j. For = 0, we have the uniform solution 


(15) 


with the density 


5 ( 1 - 7 ^) 



and the superfluid velocity Vj = hdxOj/m = const. We 
consider a countersuperflow state with a non-zero relative 
velocity 


Vr = \Vx-V2\. (17) 


Here, we will discuss the linear stability of the uniform 
countersuperflow with m = n 2 = n, as related to the 
bulk state ©. The BdG analysis gives the dispersion 
relation @,0 


If Vr > 14, the system becomes dynamically unstable 
with the non-zero imaginary part, Im w yf 0. The pa¬ 
rameter region of condition (IT^ is represented in Fig. [T] 
(b). 

The countersuperflow can be unstable even when Vr < 
14 if a collective mode causes a negative-energy fluctu¬ 
ation SE = < 0- Here, we con¬ 
sider the positive norm J2j P) > 0 without 

loss of generality. Note that dispersion (fT51) depends on 
Fg as w = oj{Vg = 0) -I- qVa- We have w < 0 with 
IF?I > Vl = ming[a;(VG = O)/^], and so the system is 
energetically unstable. This is the so-called Landau in¬ 
stability [l|. This instability is physically related to the 
motion of the center of mass of binary condensates rel¬ 
ative to the environment, such as an external potential 
or thermal excitations. The negative energy mode with 
w < 0 is spontaneously excited and amplified due to en¬ 
ergy dissipation. In this work, we do not consider the 
case where the motion of the center of mass is finite, 
since we are interested in the maximum stability of the 
Rabi-coupled countersuperflow. 

For the sake of the discussion in Sec. IIII B[ we will also 
present the stability of a countersuperflow in a finite-size 
system. What we need to show here is the maximum 
wave number qc of the unstable mode, which has Im w ^ 
0. According to Eq. we obtain 

q, = Y^Vl-V^. ( 20 ) 

The countersuperflow is stable when qc is smaller than 
1/L, where L is the system size. Note that we have 
Im w = 0 with Vg = 0 for |g| > qc and Re uj arise from 
zero at jgl = qc- This means that the system is marginally 
stable against the Landau instability for qc 1/L in the 
sense that the system possesses an energetic instability 
for any non-zero value oi Vg] Vr = Q. A characteristic 
behavior of the energetic instability of a countersuperflow 
is in the momentum change 5Jj = hq — \vj\^') ow¬ 

ing to the instability. The instability causes a relaxation 
of the relative motion between the two components, and 
the fluctuations cause a momentum change in the op¬ 
posite direction, with jgl > gc > 0; < 0. This 

characteristic behavior is revisited in Sec. IIII B I when we 
discuss the instability of the soliton. 


III. SINGLE-SOLITON STATES 


(feu - VGhqf = el + eq (^gn + 

± ^ 2 mV^{fiq + 2 gn) + (25771)2, ( 18 ) 

where I4? = (^i^iiF + rn2n2V2)/{mini + 7712712) = 
(Fl + F 2)/2 and €q = / 2 m. The countersuperflow 

is dynamically stable when 


Ffl < F = 2 


571(1-7) 


m 


(19) 


In this section, we will investigate how Rabi coupling 
affects a countersuperflow. In the presence of Rabi cou¬ 
pling, a uniform countersuperflow with = mVRx/h is 
not a stationary solution. We consider a relative velocity 
Vr = 2Trh/mL with a system size L, where the relative 
phase winds once through the system. In this case, a 
single-soliton appears as a kink of relative phase by mak¬ 
ing a density depression. Then, we discuss the stability 
of the single-soliton solution. The multi-soliton solution 
is discussed in the next section. 
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A. The single-soliton solutions 

First, we will present our numerical solutions of 
the single-soliton before performing a detailed theo¬ 
retical analysis. The solution is obtained by solving 
Eq. (1141) numerically using imaginary time propaga¬ 
tion (or the steepest descent method). The numeri¬ 
cal computations were performed from an initial state 
V'i = 7 i-(x/l+i/2) m^(jer the Neumann bound¬ 

ary condition. The system size L is taken to be large 
enough so that the system size does not affect the solu¬ 
tion. 

Figure [H shows a typical profile of the solutions. The 
relative velocity vr^x) = hdx6-/rn is localized around 
X = 0 by forming a kink of relative phase and a de¬ 
pression in the total density ni + n 2 . There is no differ¬ 
ence between the density profiles with ni = n 2 , and the 
total phase 0+ = 9i + 02 is spatially constant. 


with k ^ j. From Eq. (EU, we can write the asymptotic 
behavior of 9j (x) for x —>■ ±oo as 

9,{x) ^ (-1)^+1 , (23) 

where 

The asymptotic form of rij is derived by inserting 
Eq. (1^ into Eq. (1^ . With regards to the asymp¬ 
totic form of th e densitie s, there are three length scales: 
o-a/2, = i/ V2(l + 7 ) and . 

When (Ta/S > ^±, the asymptotic form is described as 

n,(x)-no(l-e-2|"l/'^“). (25) 

The condition (Ta/2 > reduces to 



-20 0 20 


x/^ 

FIG. 2. A typical numerical solution (black dashed curve) of 
the single-soliton with 7 = 0.4 and fl/flo = 0.04. The gray 
dashed curves show the analytical result obtained from the 
variational analysis. 


To represent the spatial profiles of the numerical results 
analytically, we performed a variational analysis, which 
is useful for simply characterizing the single-soliton solu¬ 
tions using parameters 12 and 7 . We first determined the 
asymptotic behavior of nj{x) and 9j{x) for x —00 for 
constructing a variational ansatz. Equation m is then 
reduced to 


» _ y2 d9j duj 

9x^ Uj dx dx 

h^o y/nj dx'^ 


- 9k), (21) 


dx 


'7 

no 


nk 

1 — 

no 



Ok), 


412 

12 o 


1 < 7 < 1 


312 


(26) 


which is always satisfied for stable solitons, as shown in 
Sec. IIIIBI (see also Fig. [S]). 

We assume ni(x) = n 2 (x) from our numerical results. 
If the spatial derivation of nj (x) is small so as to neglect 
the second term in the right hand side of Eq. m , one 
obtains the sine-Gordon equation 


(j 


2 

a 


2 



-b cos0_ = A, 


(27) 


where A is an integration constant, and 0+ = const, be¬ 
cause we have not considered a center of mass motion. 
The solutions of Eq. (EZl) depend on the boundary con¬ 
ditions. Under the boundary condition 0_ —>■ tt ± tt and 
dx9- ^ 0 for X —>■ ± 00 , which is equivalent to A = 1 , a 
solution is 


0_ (x) = 4 arctan , (28) 


which is called the sine-Gordon kink. This result has 
been obtained by neglecting the spatial derivation of nj 
in the limit 7—^1 [l5l |. 

By considering the asymptotic behavior in Eq. (I25L 
we constructed a variational ansatz for the density, 


nj{x) 



1 — A„sech^ — 

(Jn. 


(29) 


where A„ is the variational parameter. By inserting Eqs. 
(1251) and (1251) into the energy ([2]), and minimizing the 
energy with respect to A^, one obtains 


2012/12o 

412/12o -l- 5(1 -l- 7 ) 


(30) 


( 22 ) 


The variational ansatz proves to be a good fit to the 
numerical result in Fig. [2] 
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Here, we show the condition for applicability of our 
variational ansatz. Since the form (1281) is obtained by 
neglecting the second term in the right hand side of Eq. 
(HID, the condition is satisfied for Ay <C 1 or 

5 , , 

n;«yj(l+7)~l (31) 

under our assumption 0 < 7 < l + fl/Ho. Our approxi¬ 
mation has no strong restriction for 7 although the limit 
7 —>• 1 is assumed in Ref. (l5j| . 

To compare the analytical results with the numerical 
ones in more detail, we investigated the dependences of 
the soliton width, a, and density depression, A, on H and 
7 , which are defined by 



and 


A = 


^max ^min 


^max 


(32) 


(33) 


with max[n] = rimax and min[n] = Umin- When we use 
forms ((^ and (1^ . we have cr = CTq and A = A„. Fig¬ 
ure [3] shows the comparison between the analytic results 
(ctq, Ai,) and the numerical results ((t„, A„) of the vari¬ 
ables a and A. The two results coincide for small ft, 
consistent with condition dSB), while Ay becomes slightly 
different from A^, for large H. There is no numerical data 
in a parameter region of the plots where the single-soliton 
solution itself is unstable, which is revealed below. 




O—n/fio = 0.01 • —n/n„ = 0.05 □—n/n„ = o.07 ■ —n/n„ = o.ii 


FIG. 3. The soliton width Ca (a) and the density depression A 
(b), as functions of 7 for different Q. The numerical and ana¬ 
lytical results are displayed with marks and lines, respectively. 
The gray solid lines show the stability phase boundaries of the 
single-soliton solution, as discussed in Sec. niBi 


B. Stability of the single-soliton 

Here we investigate the stability of the single-soliton. 
According to the basic stability analysis presented in Sec. 
Ill HI it is expected that the single-soliton solution is un¬ 
stable if the maximum relative velocity ^ 2TTh/maa at 


the bottom of a density depression is large enough. We 
show here that the stability of the single-soliton states is 
explained well based on the stability analysis of counter¬ 
superflow. 

To identify the instability, we observed the dynam¬ 
ics of order parameters in the imaginary time prop¬ 
agation of our numerical simulation. The dynamics, 
which effectively show a relaxation dynamic in energy- 
dissipative systems and does not represent an actual 
time development, gives us useful information, that is, 
what kinds of mode triggers the instability. Figure 0] 
shows the instability development of the single-soliton in 
the imaginary-time propagation. The density difference 
\ni — n 2 \ starts to grow, and the momentum difference 
I Ji — J 2 I decreases at around cc = 0. Here, the local mo¬ 
mentum density Jj of the jth component is defined as 
Jj{x) = h{il;*dx4’j ~ When the density of 

a component vanishes at a given point, the kink config¬ 
uration in the relative phase is broken. After that, the 
relative velocity decays and both components flow with 
the same velocity. 



FIG. 4. A typical development from left to right of densities 
(a), phases (b), and momentum (c) in the numerical compu¬ 
tation of the imaginary time propagation. 


Since the instability starts by reducing the relative mo¬ 
mentum of the two condensates, we may expect that 
countersuperflow becomes unstable locally at the density 
depression. As a first step, we tried to apply the sta¬ 
bility criteria for a uniform countersuperflow. However, 
we found that the single-soliton solution can be stable 
even if the maximum relative velocity Umax = max[ufl] at 
the center of the density depre ssion exceeds the critical 
value Vc = 2 which is obtained simply 
by using Eq. (HU) and the minimum value nmin, of the 
density. Therefore, we need to make a correction to the 
stability analysis in order to explain the instability of the 
soliton. 

Since the instability occurs locally within the width 
^ tJo in the density depression, the finite-size-effect dis¬ 
cussed in Sec. Ill HI should be crucial for understanding 
the instability criteria of the single-soliton. In this sense. 
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the instability occurs when the soliton width (Tq is com¬ 
parable to the length 1/gcj where Qc is the upper limit 
dini). Then, we write this condition as 

qc = —■ (34) 

From Eq. (uni, and using the minimum density rimin 
and the maximum velocity r^max) the wave number qc is 
expressed as 


fn dgrimin /, ^ 

9c = ^\/<ax- (1-7- (35 

ft V m 

We use the variational results of the formula of z^max as 
a function of (Tq. Then, from Eq. (1251) . we have 

2h , , 

'^max — ■ 

m(Ta 

Therefore, the criterion of instability of the single-soliton 
is described by 


<Ja 



(1-7)- 


(37) 



FIG. 5. The stability-phase diagram of the single-soliton. The 
gray area surronnded by the black curve represents numeri¬ 
cal results of the parameter region where the single-soliton 
is stable. The gray dashed curve shows criterion dSZl) with 
ftmin = fto(l — A„). The black dashed curve shows criterion 
with numerical values of rimin. Our analysis is available 
below the gray line [see Ea. (l26ll ]. 


IV. MULTI-SOLITON STATES 


In this section we discuss multi-soliton states where 
kinks of relative phase are so close together that stable 
solutions cannot be described by the single-soliton solu¬ 
tion. The period of the soliton patterns should be 27rd, 
with which relative phase winds 


d = 


h 

mVn 


(38) 


When d is much longer than the width aa of the single- 
solitons, the pattern is described as a periodic array of 
single-soliton solutions. However, this is not true when 
the solitons are close to each other, i.e. when d ^ aa- 
Therefore, the multi-soliton states are characterized by 
the ratio 



(39) 


Under this consideration, it is straightforward to extend 
the analysis of the single-soliton solution in the previous 
section into that of the multi-soliton solution. Figure [6] 
shows the plots of Eq. ((39|) for several values of d/aa- 
Multi-soliton states appear for smaller H and larger Vr. 



n/flo 

FIG. 6. The fl-dependence of Vr for d/aa ~ 0.25, 0.5,1, 2. 


A. The multi-soliton solutions 


Figure [5] shows the stability-phase diagram of the 
single-soliton obtained from our numerical computation 
of the imaginary time propagation together with the 
plots of Eq. (021). The analytical result based on the 
variational analysis is obtained from Eq. 071) with 
rimin = zT-o(l ~ A„). The analytical result describe well 
the numerical one for smaller H, where our approxima¬ 
tion is available. The semi-analytical results are obtained 
by using numerical values of rimin in Eq. 021). The semi- 
analytical result coincides very well with the numerical 
results in Fig. [5] This result shows that the stability of 
the single-soliton states is explained quantitatively the 
stability analysis of countersuperflow by taking into ac¬ 
count the finite-size-effect. 


Numerical solutions of multi-solitons were obtained 
under the periodic boundary condition. We set the 
initial state for the imaginary propagation as i/j = 
integer number Z, where the ini¬ 
tial relative velocity was Vr = 2'Klh/mL. 

Figure 13 shows typical numerical results of the multi- 
soliton solution. When the inter-soliton spacing becomes 
of order the soliton width, d/aa ^ 1, the relative veloc¬ 
ity is no longer localized, thus making a finite relative 
velocity between the density depressions in the pattern, 
although the density profile of each depression is still sim¬ 
ilar to that of the single-soliton solution [Fig. 13 (a)]- The 
spatial profiles of the density and relative velocity be¬ 
comes similar to a sinusoidal wave for smaller d/aa [Fig. 
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FIG. 7. Numerical solutions with 7 = 0.4 and fl/flo = 0.04. 
The density difference |ni — n 2 | is zero, and the total phase 
6+ is constant (not shown). The parameters are set as 
(a) m^Vn/h = 2.51 x 10~^, where d/aa = 1.13, and (b) 
m^VR/h = 9.22 X 10~^, where d/aa = 3.07 x 10“^. 


[T^b)]. Then the density becomes substantially lower than 
the bulk value, no. 

The approximated form of the relative phase is ob¬ 
tained as a multi-soliton solution of the sine-Gordon 
equation (EH), for 2 I > 1. The solution is written as 

m 


0- (x) = 2arctan 


Vl —s? sc 


X 

- 

KUa 


where 


K = 



< 1 


(40) 


(41) 


and sc[u,m] is the Jacobi elliptic function. As a natural 
extension of the single-soliton ansatz (I29L we use the 
following form for the density anzats. 


n{x) = n„ 


1-a:: 


cn 


X 

Kaa 



(42) 


where K[77i] is the complete elliptic integral of the first 
kind, and A™ is a variational parameter of the multi- 
soliton solution. The spatial periodicity 27rc? is written in 
terms of k and Ca as 

27r(J = 2K[K]Kaa. (43) 

Here, rimax represents the maximum density between 
the density depressions. By combining Eqs. m and 
(froi) with Vi = —V 2 = Vr/2, fj,i = fi 2 = fi and 


Umin = inin[ui{] 7 ^ 0, the formula of rimax is written ap¬ 
proximately as 


1 


5(1 + 7) 

= no < 1 - 


hn 


- -:rV„ 


f m^Vn 


4(1 -f 7 ) V ^ 


(44) 


The last term in the right hand side of Eq. dm) shows 
a modification due to the non-zero relative velocity be¬ 
tween density depressions. 

The variational parameter A™ is obtained by inserting 
Eqs. dm 3(iid. (j42l) into Eq. m, and then minimizing it 
with respect to A™. Our variational ansatz dm and 
(1421) agree with the numerical results of the multi-soliton 
solutions in Fig. 0 Although our variational calculation 
can be inadequate for large H, our variational calculation 
is available in a wide range of parameters in Fig. [SJ since 
the multi-soliton states are stable for relatively smaller 
values of H as is disscused below. 


B. Stability of the multi-soliton 

Stability analysis of the single-soliton states were per¬ 
formed by considering the stability of countersuperflow 
around a density depression. For stability analysis of 
multi-soliton states, we have to consider an additional 
possibility that instability occurs in high-density regions 
because of nonzero relative velocity there. 

To obtain a stability criterion for the former possibil¬ 
ity, we discuss the stability of multi-soliton states in a 
manner similar to the stability analysis for single-soliton 
states. From Fa. (HOI) , the width of a density depression in 
the multi-soliton solution is represented by a length kctq . 
Instability can occur when \/noa equals the critical wave 
number (1^01) . The phase boundary of the instability ow¬ 
ing to the density depressions is written in a similar form 
to Eq. (I37l) as 


— = J^(l-7). (45) 

KOa V 3no 

This criterion reduces to the stability criterion (1571) in 
the single-soliton limit djoa —>■ 00 with Kaa a- In 
the opposite limit d/aa —>■ 0 with d fixed, density de¬ 
pressions are negligibly small with H/Ho —>■ 0. Then, the 
criterion (|45ll must be reduced to the criterion of uniform 
countersuperflow obtained from Eq. (I19p with I 4 = Umax 
and n = rimin- However, the criterion (1451) in the limit 
does not give a correct results owing to the asymptotic 
behavior naa —>■ 2d for d/aa 0, with which the finite- 
size-effect remains even for H = 0 without density de¬ 
pressions. 

This inconsistency is recovered by the additional pos¬ 
sibility that instability occurs in high-density regions, 
where the density takes the maximum value, n = rimax- 



















Supposing that instability occurs when the relative veloc¬ 
ity Umin in high-density regions exceeds the critical value 
(fT^ with n = Umax, the stability criterion for small fl/flg 
is written as 

— = ./^(l-7). (46) 

K(Ta V 2no 

This criterion reduces to that of uniform countersuper¬ 
flow consistently in the uniform countersuperflow limit, 
d/aa —>■ 0 and il/ilo 0. 
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FIG. 8. (a) The stability-phase diagram of the Rabi-coupled 
countersuperflow. The gray solid curve shows the stability- 
phase boundary in Fig. [S] The gray dashed curve represents 
the stability criterion m of uniform countersuperflow, (b) 
The two dimensional plot of (a). The black dashed curves 
show the criterion (|45l) with numerical data of rimin. The 
gray dashed curves show the criterion dSli. 


Stability-phase diagram of our system is summarized 
in Fig. [H Numerical plots are consistent with the semi- 
analytical estimations based on Eqs. (HSl) and (H51) . The 
analytical curve of Eq. (gSl) crosses the semi-analytical 
curve of Eq. (I45p in the left side of the phase diagram 
in Fig. [8](b). This structure is also consistent with the 
numerical plots. These results show that the stability of 


this system is described totally by the stability analysis 
of countersuperflow. 

V. SUMMARY AND DISCUSSION 

We studied the Rabi-coupled countersuperflow of bi¬ 
nary Bose-Einstein condensates in quasi-one-dimension. 
The variational formulas provide a good description of 
the stationary states of the single-soliton [Eqs. (l28l) and 
d^ ]. and the multi-soliton solutions [Eqs. (HUl) and (H^ ] 
for small Rabi frequencies. By taking into account the 
finite-size-effect due to density depressions in the soli- 
ton solutions, the stability analysis of countersuperflow 
is applicable when explaining the stability-phase diagram 
of the Rabi-coupled countersuperflow, Fig. [S] for single- 
soliton states, and Fig. [8] for multi-soliton states. 

These solitons will be observed in experiments as a 
density contrast in atomic clouds. Therefore, the param¬ 
eter dependence of the density depression in Fig. [5] is an 
important benchmark for observing solitons experimen¬ 
tally. For a density contrast of higher than 5%, A > 0.05, 
7 = 512/5 must be smaller than ^ 0.9. This condition 
will be achieved in future experiments by utilizing Fes- 
hbach resonance [ 2 ^. Even for 7 ^ 1 , we may expect 
a unique behavior of this system; e.g. instability de¬ 
velopment can be distinct from those of countersuper¬ 
flow instability without Rabi-coupling [E H' B and Rabi- 
coupling-dominant pattern formation [l0l[ . Studies of in¬ 
stability development provides an interesting framework 
for future investigations. 
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